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ABSTRACT 

In this paper we investigate the gravothermal instability of spherical stellar systems 
endowed with a radially anisotropic velocity distribution. We focus our attention on the 
effects of anisotropy on the conditions for the onset of the instability and in particular 
we study the dependence of the spatial structure of critical models on the amount of 
anisotropy present in a system. The investigation has been carried out by the method 
of linear series which has already been used in the past to study the gravothermal 
instability of isotropic systems. 

We consider models described by King, Wilson and WooUey-Dickens distribution 
functions. In the case of King and WooUey-Dickens models, our results show that, for 
quite a wide range of amount of anisotropy in the system, the critical value of the 
concentration of the system (defined as the ratio of the tidal to the King core radius 
of the system) is approximately constant and equal to the corresponding value for 
isotropic systems. Only for very anisotropic systems the critical value of the concen- 
tration starts to change and it decreases significantly as the anisotropy increases and 
penetrates the inner parts of the system. For Wilson models the decrease of the con- 
centration of critical models is preceded by an intermediate regime in which critical 
concentration increases, it reaches a maximum and then it starts to decrease. The 
critical value of the central potential always decreases as the anisotropy increases. 

Key words: globular clusters: general - stellar dynamics 



1 INTRODUCTION 

It is well known that the phenomenon of gravothermal in- 
stability is strictly related to the property of self-gravitating 
systems of having a negative specific heat. In 1962 Antonov 
showed that there is no state characterized by a local maxi- 
mum in the entropy for stellar systems with fixed mass and 
energy, contained in a spherical box with radius greater than 
the critical value rcrit = -0.335{GM'^ /E) where M is the 
total mass of the system and E its total energy. Lynden-Bell 
and Wood (1968) interpreted this result in their pioneering 
work on the stability of isotropic isothermal stellar systems 
truncated in radius, establishing the conditions for the onset 
of the phenomenon, called then gravothermal catastrophe. 
The investigation was carried out by means of the method 
of linear series formulated for the first time by Poincare 
(1885); this technique is connected with the study of global 
stability of the equilibria of a generic system. 

In the case of isothermal configurations, the equilib- 
rium states in the presence of encounters are determined 



by finding the conditions for a maximum of the total en- 
tropy S, keeping fixed the total energy and mass, E and M, 
of the system. The phase-space distribution function that 
describes a system under these requirements turns out to be 
the Maxwell-Boltzmann function, 

/(£)=Ae-''^ (1) 
where 

e=^v^ + ^r) (2) 

is the energy per unit mass of single stars, (3 is the inverse 
of the "temperature", ^ is a normalization constant and 
$(r) is the self-consistent gravitational potential obtained 
by solving the Poisson equation. Lynden-Bell and Wood fo- 
cussed their attention on the study of the stability of a sys- 
tem against perturbations keeping constant, besides the to- 
tal mass and energy, the radius of the system also. The criti- 
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cal value of the concentration of the system, as measured by 
the ratio of the central density, pc, to the average density in- 
side the boundary of the system, pb, is equal to pc/pb = 709 
which coincides with the result found by Antonov. 

An analogous investigation was later carried out by 
Katz (1978, 1980) for more realistic models of stellar clus- 
ters. Katz studied the gravothermal instability of King 
(1965), WooUey-Dickens (1961) and Wilson (1975) models, 
which are all characterized by distribution functions with a 
cut-off in the energy space. In this case, the study of the 
stability was realized keeping fixed the total mass, the to- 
tal energy and the multiplicative normalization parameter A 
which reflects the indetermination in the zero of the chemical 
potential (see the appendix of Katz (1980) for the discussion 
on the motivation for this choice). In these cases, the onset 
of gravothermal instability occurs for the following values of 
the central dimensionless relative potential, 'i'c = P'^c, 



1) \['c ~ 7.65, Woolley — Dickens models, 

2) = 7.40, King models, 

3) ^'c = 6.67, Wilson models. 



(3) 



Taking anisotropy in the velocity distribution into ac- 
count is an essential step in order to improve the reliabil- 
ity of energy-truncated models to better describe the struc- 
ture and evolution of clusters. The important role played by 
anisotropy has been recognized since the dawn of the study 
of globular clusters (Woolley & Robertson, 1956). The main 
goal of the present paper is that of extending the study 
of gravothermal instability to spherical systems endowed 
with an anisotropic velocity distribution and determining 
the dependence of the critical value of VPc on the amount 
of anisotropy present in the systems. We point out that our 
analysis is meant to provide only a general indication of the 
role of radial anisotropy on the onset of gravothermal in- 
stability, providing a simple and straightforward estimate 
of its effect as measured by the critical value of the cen- 
tral potential and of the concentration of the system (given 
by the ratio of the total radius of the system to the King 
core radius). A much deeper insight in the collisional evo- 
lution of anisotropic systems can be obtained only by other 
methods of investigation such as TV-body simulations, in- 
tegration of the Fokker-Planck equation and gaseous mod- 
els (see e.g. Takahashi 1995, 1996, Takahashi, Lee & In- 
agaki 1997, Spurzem & Takahashi 1995, Giersz & Spurzem 
1994, Bettwieser & Spurzem 1986, Giersz & Heggie 1994a,b, 
1996, 1997, Louis & Spurzem 1991, Spurzem 1991, Spurzem 
& Aarseth 1996, for some interesting and detailed investi- 
gations on the evolution of anisotropic stellar systems). In 
Sects. 3.3 and 4, we compare the global indications resulting 
from the present study with the specific findings of some of 
the most recent works cited above. 

The scheme of the paper is the following; in Section 2 
we summarize the main properties of these models, in Sec- 
tion 3 we study their stability by means of the method of 
linear series discussing the results and, finally, in Section 4 
we summarize the main conclusions. 



2 ANISOTROPIC MODELS 

The relevance of an anisotropic model for the description of 
real stellar systems arises when studying the first periods 



of life of stellar systems: it is possible to demonstrate that, 
at the end of the process of the violent relaxation (Lynden- 
Bell, 1967), stellar systems are anisotropic because they still 
have stars mainly distributed on radial orbits, and, although 
they tend to releix and become isotropic on a time scale 
which is of the order of the binary relaxation time (see, 
e.g., Binney & Tremaine, 1987), if the conditions for the 
onset of the gravothermal instability take place before their 
total isotropization, anisotropy cannot be left out of con- 
sideration. It must be remarked that, in evolving systems 
like globular clusters, anisotropy is not simply removed by 
relaxation. Rather, it is even increased in the halo, due to 
the high energy stars coming from the strong relaxing core 
(Takahashi 1996) and it is also expected to gradually pene- 
trate into the inner regions as the core collapses. The situa- 
tion is different for tidally truncated systems in which stars 
can escape at a significant rate from the cluster: in this case 
it has been shown (Giersz & Heggie 1997, Takahashi et al. 
1997) that anisotropy production can be partially or com- 
pletely balanced by the evaporation of stars which would 
occur preferentially for stars on radial orbits. 

In any case, if primordial anisotropy changes the con- 
ditions for the onset of gravothermal instability, this can 
also have important consequences on the subsequent evolu- 
tion of anisotropy itself. In the present Section we describe 
the properties of the anisotropic variants of the energy- 
truncated models mentioned above, which, in the subsequent 
one, will be submitted to the analysis of gravothermal sta- 
bility. 



2.1 Description of anisotropic spherical models 

Following arguments similar to those of Tremaine (1986) in 
the discussion of the process of violent relaxation, the distri- 
bution functions introduced to describe anisotropic systems 
are obtained by multiplying by a factor of the form 



7^ 



(4) 



the isotropic King, WooUey-Dickens and Wilson distribution 
functions, which therefore can be written in the form 



1) fwD{e,L')^Ae-^ e 



(e-*-l) 



(5) 



3) fw{e,L^) = A e ^ + 



The above definitions hold for e < 0, while it is assumed that 
/(e) = for e > 0. In each distribution function, the param- 
eter /3 has been substituted with the inverse of the isother- 
mal velocity scale . The dependence of the functions on 
the square of the angular momentum L ensures the spherical 
symmetry of the global properties of the models. In (^), ra 
is the anisotropy radius that represents the value of the ra- 
dius beyond which the orbits start to be mainly radial. The 
distribution function /i<-(e, L^) is usually known as Michie 
model (Binney & Tremaine, 1987), whereas fw{e, L^) is sim- 
ply the spherically symmetric version of the generic function 
introduced by Wilson (1975) himself. 
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In order to perform the integration of these models, we 
define the following dimensionless quantities: 



X 



r 

rK ' 
so that 



V 

a 



In 



e 1 -2 5. 



3a 



$ = — . 



(6) 



(7) 



(8) 



is the Kmg radius and gives the length scale. Calling rt the 
tidal radius of the system, so that Xt = rt/rK, the self- 
consistent relative potential 



*(a::) = -$(a;) + <I>(a;t) 



(9) 



and, from it, the physical properties of the models, is ob- 
tained by solving the Poisson equation which can be writ- 
ten in terms of the dimensionless quantities introduced in 
(|) and (0) as: 



1 d 

x^ dx 



dx 



d^* 2 d* 



da;2 X 



dx 



^p{'i>{x),x)) 



-9 



(10) 



where the last equality comes from the substitution of tk 
with its explicit expression. The related Cauchy's problem is 
then solved with initial conditions *(0) = *c and *'(0) = 0. 

The degree of anisotropy of a given model is 
parametrized by the dimensionless anisotropy radius 



7 



rK ' 



(11) 



We have carried out the integration of anisotropic WooUey- 
Dickens, King and Wilson models by selecting different val- 
ues of the anisotropy parameter 7 and of the central poten- 
tial ^c. 7 has been chosen to range from almost isotropic 
configurations (7 = 100) up to highly anisotropic models 
(7 = 0.3). Models with higher values of ^I'c tend all to the sin- 
gular isothermal sphere; however, highly concentrated mod- 
els, pose a certain constraint on the range allowed to the 
anisotropy parameter. On the other hand, following argu- 
ments analogous to those in Merritt et al. (1989), it can be 
shown that singular models with a completely radial distri- 
bution display unphysical behaviour. 



2.2 Infinite models 

An interesting feature of these models is the existence of def- 
inite pairs of values (7, vE'c) that identify unbounded models, 
i.e. models with an infinite radius. In fact, there is a criti- 
cal value of the central potential, depending on the degree 
of anisotropy characterising each model, beyond which we 
have rt —> 00 for every 7; this happens even though the dis- 
tribution functions used for the description of these models 
are energy-truncated and the corresponding isotropic config- 
urations are always bounded. Such remarkable peculiarity is 
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Figure 1. 733 versus the dimensionless central potential 'I'c for 
the first unbounded Woolley-Dickens (a), King (b) and Wilson 
model (c). 



common to all anisotropic models, regardless of the form of 
the dependence on the energy of the functions {js]). 

Following Stiavelli & Berlin (1985) we consider the pa- 
rameter 

1 



7BS 



(12) 



4TTGAari ' 

This parameter is related to our anisotropy parameter 
through the following relation: 

Pc 



7BS — 



(13) 



In Fig. 1 we have plotted 7bs versus 'I'c for the 
anisotropic Woolley-Dickens, King and Wilson models hav- 
ing an infinite radius. These models are singled out identi- 
fying, along the sequence at given ^'c parametrized by 7bs, 
the first system for which the relative potential goes as 1/x 
at large radius. 

It is evident that the behaviour is the same of that 
shown in the paper by Stiavelli & Berlin (1985) and thus 
it seems that this characteristic is common to all the spher- 
ical anisotropic models having a dependence on the angular 
momentum of the kind e , regardless to the form of the 
energy distribution. 

Note that not all these unbounded models are rele- 
vant to model real physical systems; in fact, while the total 
dimensionless mass of infinite Wilson and Bertin-Stiavelli 
models is finite, it diverges for anisotropic Woolley-Dickens 
and King models, making the latter models useful only for 
those values of and 733 which lead to a finite-size system. 



GRAVOTHERMAL INSTABILITY OF 
SPHERICAL ANISOTROPIC MODELS 
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3.1 Linear series of equilibrium 

The investigation of the gravothermal stability of the mod- 
els discussed in Section 2 has been carried out by means of 
the method of linear series. This approach, introduced by 
Poincare in 1885, is based on the study of global stability 
of the equilibrium states of a generic system. In many sta- 
tionary problems, the description of a system depends on 
some parameter other than the generalized coordinates qi 
{i = 1, . . . , n) in terms of which we describe the configura- 
tion of the system. If we call fi this parameter and U the 
potential energy of the system, the condition for the ex- 
istence of an equilibrium is that the tangent plane to the 
n- dimensional surface 

U{fi, Qi) = const (14) 

in the equilibrium point is perpendicular to the fi axis. For 
every fixed value of ^ in appropriate ranges, we obtain an 
equilibrium configuration. Therefore, if we call this con- 
figuration, we have 

g° = g°(M). (15) 

These points are called level points and, by varying the value 
of the parameter /i, they form a curve called linear series 
(Poincare, 1885). To verify the stability of these equilibria, 
we restrict our attention to a plane fi — const. In statics, 
the condition for an equilibrium to be stable is that the 
potential energy U in that point has a minimum. In our 
case, this statement means that a point on the linear series 
( p^ represents a stable configuration only if all the different 
vertical sections of the surface U = constant across that 
point are turned towards the direction of decreasing U. 

It is possible to prove (see, e.g., Jeans, 1928) that ev- 
ery change in the stability of a system corresponds to an 
extremum of the linear series. Accordingly, the investigation 
of the stability of a self-gravitating system with respect to 
the gravothermal catastrophe, is accomplished locating the 
first critical point on a suitably parametrized linear series. 
This is a sequence of equilibrium systems with appropriate 
constraints insuring that the perturbation implicit in the 
stability analysis, keeps fixed the relevant physical quantities 
of the system. The critical value of the parameter identifies 
the transition from stability to instability. 

3.2 Choice of the constraints 

In the problem we are investigating, since the distribution 
functions ^ depend on the three parameters A, a, 7 (or 
Ta) and on the central potential ^Jc, we will have to fix, 
in agreement with the theory, three quantities belonging to 
the particular model under consideration and let a fourth 
quantity vary as a function of a parameter that is able to 
describe such model in a unique way. Besides keeping fixed 
the total mass of the system and the parameter A (see Katz, 
1980, for a discussion about this constraint), the parameter 7 
has been held constant and we will investigate the stability 
by means of the curves of the total energy of the system 
versus the central dimensionless potential. Our results will 
thus be relevant for perturbations in which A, M, E and 7 
remain constant. 

As for the first three quantities we thus adopt the same 
constraints used in the study of stability of isotropic sys- 
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Figure 2. E vs 'I'c for WooUey-Dickens models. 



tems by Katz (1980). For what concerns the fourth quantity, 
it is worth noting that the anisotropy parameter 7 is only 
one of the possible choices, but other quantities related to 
the anisotropy of the configurations could be held constant 
along the equilibrium sequences. We could not find any con- 
vincing physical criterion in favor of one particular choice 
and the parameter 7 chosen has the advantage of allowing 
a straghtforward implementation of the calculation. 

While it would be clearly desirable that the arguments 
supporting the choices of the constraints (in particular the 
choice of keeping A fixed) were substantiated by a more rig- 
orous theory on the statistical mechanics of the gravitating 
systems, there are, however, several independent evidences 
supporting the validity of the basis of the present approach. 
In fact, in the case of isotropic King models, the results of 
some Fokker-Planck simulations (which are free of any as- 
sumption of the kind made in our investigation) by Wiyanto 
et al. (1985) have confirmed with good accuracy the insta- 
bility limit found by Katz (1980) and the Fokker-Planck 
simulations by Cohn (1980) have shown that, for Plummer 
models also, the value of the central dimensionless poten- 
tial at the onset of gravothermal instability is similar to 
that found by Katz. Moreover, it is worth emphasizing that 
the choice of the third parameter to be kept constant along 
the series is likely not to be so important for the conclu- 
sions one can draw by means of this analysis. In fact, for 
example, if in the analysis of the stability of isotropic King 
models, instead of keeping fixed A, another quantity like the 
tidal radius were kept fixed, the critical value of the central 
potential would not change significantly, being in this case 
*c =i 8 (see, e.g., the plot of E/{GM'^/rt) versus Wo = *c 
in Chernoff et al. (1986) from which this critical value can be 
obtained). Analogous results have been obtained in a more 
recent analysis by Lagoute & Longaretti (1996) who, with 
the same approach adopted in the present work and the same 
choice of the constraints, have investigated the gravothermal 
instability of rotating King models. 
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Figure 3. E vs 'I'c for King models. 







1 1 1 1 


7 0.3 ; 
















7=30/ / 


- V 


= 3 






i 

7=10| 






1 

1 


y=100 












1 , , , 1 - 



2 4 6 8 10 



Figure 4. E vs for Wilson models. 
3.3 Transition to instability 

Figg. 2-4 show the plots of E versus for different values 
of 7, for WooUey-Dickens, King end Wilson models respec- 
tively. The anisotropy parameter ranges from 100, which 
corresponds to nearly isotropic models, to 0.3, which cor- 
responds to models with a strong anisotropy. The minima 
of these curves show the configurations in correspondence 
of which we have the transition between stable (to the left 
of the minimum) and unstable models (to the right of the 
minimum). 

For large values of 7, only the very outer regions of 
the systems are affected by the presence of anisotropy and 
thus the equilibria are very similar to those of the corre- 
sponding isotropic systems. As expected, the critical values 



of tends to the critical value for isotropic systems listed 
in (|^. As 7 decreases, the anisotropy extends deeper and 
deeper into the core, changing significantly the properties of 
the equilibria so that the critical value of the central poten- 
tial is quite different from that of isotropic systems. This is 
clearly shown in Fig. 5 in which we have plotted the critical 
values of '^c versus 7 for the models investigated. The in- 
teresting point to note is that the main effect of anisotropy 
is that of causing the system to become unstable at lower 
values of 'I'c. 

However, we have to remark that systems with different 
values of 7 and the same value of do not have the same 
structure. This means that, although Fig. 5 gives an inter- 
esting result on the dependence of the critical value of the 
central potential on the amount of anisotropy in the system, 
it does not provide quantitative information on the issue of 
the change of the spatial structure of the critical model as 
a function of the parameter 7. In order to investigate this 
point, we have focussed our attention on the concentration 
parameter 

c = log , (16) 

\rK/ 

and in Fig. 6 we have plotted the concentrations of the crit- 
ical models, Ccrit, versus the parameter 7. 

We see that, for King and WooUey-Dickens models, Ccrit 
is approximately constant over a wide range of values of 7 
and it starts decreasing only for log 7 < 0.3, that is for 
systems in which the very inner regions also are signifi- 
cantly affected by anisotropy. This means that, for King and 
WooUey-Dickens models, the structure of the critical models 
is approximately independent on 7 and the conditions for the 
onset of gravothermal collapse are not significantly affected 
by the amount of anisotropy in the system; only for fam- 
ily of systems characterized by extreme anisotropy, deeply 
penetrating into the inner parts, the structure of the critical 
model changes and the instability sets in for lower values of 
c. 

The situation is different for Wilson models: for these 
systems Ccrit has a maximum for log 7 ~ 0.5 and then de- 
creases for log 7 < 0.5; in this case also, for most values of 
7 investigated, the spatial structure of the critical equilibria 
does not vary much with the anisotropy. 

While an exact comparison with the results of numerical 
simulations of anisotropic systems (see the papers quoted in 
the Introduction) is not easy, we point out that our results 
are in good qualitative agreement with the results obtained 
by integration of the Fokker-Planck equation by Takahashi 
(1996) and Takahashi et al. (1997). As shown in these works, 
in fact, the main difference in the evolution towards core 
collapse of isotropic and anisotropic systems is in the evolu- 
tionary timescale, the former evolving faster than the latter. 
There is, instead, no significant difference in the structures 
of isotropic and anisotropic systems, when compared at the 
same evolutionary phase (see Fig. 6 in Takahashi (1995) and 
Fig. 2b in Takahashi et al. (1997)). As for the work of Taka- 
hashi et al. (1997), it is important to note that the max- 
imum amount of anisotropy is reached at the moment of 
maximum contraction at the end of the process of core col- 
lapse and the best-fitting anisotropic King model at that 
phase has 7 = 300. We point out that our critical mod- 
els refer to the structure of the systems at an earlier phase 
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log 7 



Figure 5. Critical value for the onset of gravothermal instability 
of the dimcnsionless central potential versus 7 for King (dashed 
line), WooUcy-Dickens (solid line) and Wilson (dot-dashed line) 
models. For any value of 7, models with \E'c larger than the critical 
value are unstable against gravothermal collapse. 




log 7 



Figure 6. Critical value (for the onset of gravothermal instabil- 
ity) of the concentration parameter, c, versus 7 for King (dashed 
line), WooUey-Dickens (solid line) and Wilson (dot-dashed line) 
models. For any value of 7, models with c larger than the critical 
value are unstable against gravothermal collapse. 



than that of maximum contraction (our critical models are 
meant to describe the system just at the moment of onset of 
gravothermal instability), when the value of 7 in the simu- 
lations of Takahashi ct al. must be even larger. This moans 
that the onset of gravothermal instability in the simulations 
of Takahashi et al. occurs well in the regime of high values 
of 7, where our analysis shows that no difference is expected 



between the spatial structure of the critical models for the 
isotropic and anisotropic case. 

Additional simulations of systems with a stronger 
anisotropy at the moment of the onset of core collapse would 
be extremely interesting in order to test our conclusions 
about the decrease of the critical value of c for small val- 
ues of 7. 



4 CONCLUSIONS 

Models based on energy truncated distribution functions are 
still a very useful tool to explore some of the properties 
exhibited by stellar systems when they are in an enviroment 
that can heavily influence their characteristics: this is the 
case of globular clusters, but can be of relevance even in the 
case of elliptical galaxies in rich clusters. 

In this paper, after having described their main struc- 
tural properties, we have investigated, by the method of 
linear series, the stability against gravothermal collapse of 
models described by King, Wilson and WooUey-Dickens 
anisotropic distribution functions. The amount of anisotropy 
of a system has been quantified simply by means of the ra- 
tio, 7 = ra/rK, of the anisotropy radius ra, which is the 
radius beyond which radial orbits dominate, and the King 
core radius rx. 

We have shown that, as 7 decreases (more anisotropic 
systems), the central potential of the critical model de- 
creases. The spatial structure of the critical model, globally 
described by the concentration parameter c = log(rt/rK), 
remains unchanged for a wide range of values of 7 (this 
is true in particular for King and WooUey-Dickens models, 
while the range is smaller for Wilson models); in fact the 
concentration of the critical models as a function of 7 is 
approximately constant and equal to the typical value for 
isotropic systems down to log 7 > 0.3 (log 7 > 1 for Wilson 
models). Only for very anisotropic systems, with log 7 < 0.3, 
the concentration of the critical models decreases as 7 de- 
creases and the anisotropy reaches the very central regions 
of the systems. 

We have shown that our conclusion for high values of 7 
is consistent with the results of a recent investigation car- 
ried out by means of the integration of the Fokker-Planck 
equation for anisotropic systems by Takahashi et al. (1997), 
showing that no significant difference exist between the spa- 
tial structure of isotropic and anisotropic systems when they 
are compared at the same evolutionary pheise. 

We point out again that numerical simulations designed 
to follow the evolution of systems characterized by a stronger 
anisotropy at the moment of the onset of core collapse would 
be suitable to verify our conclusions concerning the decrease 
of the critical concentration for smaller values of 7. 
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